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ABSTRACT 


In the future, astronauts will be sent into space for longer durations of time compared to 
previous missions. The increased risk of exposure to dangerous radiation, such as Galactic 
Cosmic Rays and Solar Particle Events, is of great concern. Consequently, steps must be taken to 
ensure astronaut safety by providing adequate shielding. In order to better determine and verify 
shielding requirements, an accurate and efficient radiation transport code based on a fully three 
dimensional radiation transport model using the Green's function technique is being developed. 

INTRODUCTION 

The exposure of space travelers to radiation is determined in part by the transport properties 
of the radiation throughout the spacecraft, such as its onboard systems and the bodies of the 
individuals themselves. Meeting the challenge of future space programs will therefore require 
accurate and efficient methods for performing radiation transport calculations to determine and 
verify shielding requirements. According to a recent National Research Council Report [1], 
predictions derived from radiation transport calculations need to be tested using a common code 
for laboratory and space measurements that has been validated with accelerator results. Studies 
by Wilson et al. [2,3] have identified Green's function techniques as the likely means of 
generating efficient high charge and energy (HZE) shielding codes that are suitable for space 
engineering and are capable of being validated in laboratory experiments. In particular, the 
Green's function approach has the advantage of high computational efficiency compared to 
Monte Carlo codes in both laboratory and space simulations. Consequently, a laboratory code 
designed to simulate the transport of heavy ions through a single layer of material was developed 
[2,4], It was based on a Green's function model as a perturbation series with non-perturbative 
corrections. The code was validated for single layer targets and then extended to handle multi- 



layer targets [5,6,7]. This early code used a scale factor to equate range-energy relations of one 
material thickness into an equivalent amount of another material, and proceeded to perform the 
transport calculations in the new material [8]. While this method has proven to be acceptable 
using low-resolution detectors [6,9], it is not the most accurate reflection of different material 
properties and is not well suited for high-resolution measurements. Range and energy straggling, 
multiple Coulomb scattering, and energy downshift and dispersion associated with nuclear 
events were lacking from prior solutions. In recent publications [10,1 1], it has been shown how 
these effects can be incorporated into the multiple fragmentation perturbation series leading to 
the development of a new Green's function code GRNTRN (a GReeN's function code for ion 
beam TRaNsport). GRNTRN has proven to be accurate in modeling ion beams for a single layer 
of material [9,12], and has been extended to handle multiple layers [13]. Unlike the earlier 
Green's function code, the multi-layer GRNTRN code does not make use of range scaling, but 
instead transports ions through the target layer by layer. It is, however, a purely one dimensional 
code and does not take the variation of particle flux with angle into account. In order to remove 
this deficiency, it will be necessary to develop a fully three dimensional GRNTRN code. In this 
work, the progress made toward developing a fully three dimensional version of GRNTRN is 
reported. 

TRANSPORT THEORY 

Consideration is given to the transport of high charge and energy ions through a three- 
dimensional convex region V, which is bounded by a smooth surface dV, and is filled with a 
target material. According to Wilson [14], the transport process is governed by the continuous 
slowing down, linear Boltzmann equation 
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with a boundary condition of the form 

<t>j ( x b ,QE) = Fj ( x b ,QE), 


( 2 ) 


where x b is a point on the boundary, and N is the number of ions being transported. In equation 
(2), and throughout the rest of this paper, it is assumed that all boundary particles are directed 
inward towards the volume. In equations (1) and (2), F.(x b ,£l,E ) is a prescribed function, 


< p.(x,Sl,E ) is the flux of j-type ions with atomic mass Aj at position x b moving in the direction 
£2 with kinetic energy E, Oj{E) is the macroscopic absorption cross section, S-(E) is the ion 


stopping power, and a jk (El,El\E,E') is the double differential production cross section for 

interactions in which j-type ions with energy E and direction £2 are produced by k-type ions with 
energy E' and direction £2'. Throughout this paper, the units of energy will be expressed in 

MeV/amu, and units of depth in g/cm . The latter is obtained by multiplying the linear depth in 

• • t # 1 % m 
cm by the volumetric density of the material in g/cm . Unless otherwise stated, all energies are 

kinetic. 

The method of characteristics can be used to show that the Boltzmann equation and its 
boundary condition are equivalent to the transport integral equation [14] 
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where p = x • £1 is the projection of the position vector onto the direction vector, x n = x — p£l 

is an intermediate vector, and x' = x n + p'£2 is the point where the ray through x in the 

direction £2 enters V. The term R/E) is the range of a j-type ion with energy E, and the residual 
energies 


E j = R- l [R,(E) + p - p'}. 
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are obtained through the usual range-energy relationship. Lastly, the nuclear survival 
probability, Pj(E), is defined as 


Pj(E) = exp 
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By introducing the field vector <D(x,Q ,E), whose components are the j-type particle fluxes 
<pj(x,£l,E), and following the procedure described in [10], the transport integral equation may 

be expressed in the operator form as 

<D = G°-F + Q- L- E- <D. (7) 


Since equation (7) is a Volterra type integral equation, it admits the Neumann series solution [10] 
<D = [G° + (Q • L • E) • G° + (Q • L • S) 2 • G° + (Q • L • S) 3 • G° + ...] • F 

( 8 ) 

= [G° + G 1 + G 2 + G 3 +...]• F, 


where F is the boundary flux vector whose components are the j-type boundary fluxes Fa . The 

term G" for n = 0, 1, 2, 3,... appearing in this expansion is called the n th order Green's operator 
and is associated with the n generation of fragments produced. The above formalism lends to 
the following interpretation of the Neumann series. In the first term, G° • F, the operator G° 
propagates primary ions from the boundary to the point x with attenuation processes. In the 



second term, Q • L • S • G° • F , the operator G° propagates primary ions from the boundary to an 

interior point x " where a nuclear event takes place. The term E ■ G° • F is the production 
density of first generation secondaries at position x " . These ions are transported from x " to x 
by the linear transport operator L. Lastly, the quadrature operator, Q, sums up all of the first 
generation fragments that are produced at interior points and transported to the point x . The 
remaining terms of the Neumann series may be interpreted in a similar way and, from the series, 
the n order Green's operator can be computed recursively for n > 1 as 

G” = (Q • L • 3) • G” -1 . (9) 

When the boundary condition (2) takes the special form 

- n • qmb - E 0 )J(x b - jo), (io) 

where 6 is the surface delta function on the boundary, the solution of equation (1) is called the 
Green's function and is denoted by the symbol G • (x,x Q ,£l,£l 0 ,E,E 0 ) . Once the Green's 

function is known, the solution for an arbitrary boundary condition, as in equation (2), can be 
obtained from the formula 

Q E) = [G • V]j(x,£l,E) = [(G° + G 1 + G 2 + G 3 H ) • V]-(x,£l,E) 
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The summation is taken over k >j (instead of k > j) to account for the primary ion spectrum. 

ZERO ORDER GREEN’S FUNCTION 

The zero order Green's function is the first term in the Neumann series (8) with the unit 
boundary condition (10). When taking energy straggling in to account, as described in [10] and 


[1 1], it can be shown that the zero order Green's function takes the form 
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where P m (E) is the survival probability for m-type ions of energy E, and the residual energy is 

E m = R m - 1 lRJ E ) + P-Pl 03 ) 

The term 

EJp - p\E 0 ) = - (P - P% (W) 

is the mean energy at depth {p — p') of an m-type ion that entered the transport material with 
kinetic energy E 0 , and s m (p — p\E 0 ) is the corresponding energy straggling width. 

When the boundary condition takes a more general form, as in equation (2), the primary flux 
is obtained by integrating the zero order Green's function over all energies and directions within 
the volume 

^°(x,£J,B) = [G°-F],(x,n,E) 

/ ■ r noo 0 (15) 
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which, in general, has to be evaluated numerically. In the accelerator beam model described 
below, equation (15) can be approximated analytically and a closed form expression may be 
obtained for the zero order primary flux. The result obtained in this case is called the broad zero 
order Green's function. 

Since ion beam experiments play an important role in analyzing the shielding requirements 
against space radiations, modeling the propagation of accelerator beams through potential 
shielding materials is of interest. A simple model can be constructed by assuming that the 
accelerator beam consists of mono-energetic m-type particles of energy Eq, that move in the 



direction £2 q and enter the material at the point Xg on its boundary. In this case, the boundary 


condition is given by equation (10), and the solution is the Green's function. In practice, 
however, the accelerator beam is neither mono-energetic nor unidirectional, and it enters the 
material at several points that are clustered about the mean value x 0 . A more realistic model can 
therefore be obtained by assuming that the beam has Gaussian profiles in both angle and energy, 
and that it enters the material at points that are distributed in a Gaussian manner about the mean 
point of entry, x 0 . In order to accomplish this, we will assume that the boundary is defined by 
the single-valued, continuously differentiable parametric equations dV = {x : x = x(m,u)}, 
with the bounds for u and v given by u s <u<Uj,v s <v<Vj. In this case, the element of 
surface area is given by dS = \d u x x d v x\dudv . Further, at a non-singular point, 

x Q = x 0 (u 0 ,v 0 ), of the parametric coordinate system, the surface delta function is denoted by 8 
and is given by 

8(x- x 0 ) = |3„x 0 x d v x 0 1 _1 8(u - u 0 )8(v - v 0 ) . (16) 


The boundary condition (2) may then be assumed to take the Gaussian form 
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where x b = x b (u b ,v b ), H[x] is the Heaviside Function, and s x , sq, and s E are the spreads in 
space, angle and energy, respectively. The normalization constants K n and K x are given by 
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where erf refers to the error function. It should be observed that in the limit as s x , sq, s e — > 0 , 

the boundary condition (17) reduces to the Green's function boundary condition (10). 

Equations (12) and (17) may now be substituted into equation (15) and the resulting integrals 
approximated by the mean value theorem and saddle point techniques discussed in [10]. The 
primary flux, which in this case is called the broad zero order Green's function, 

G b jm (x, x 0 ,n, Sl 0 ,E,E 0 ), is then given by the expression 
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where the residual energy is 
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and the broad energy spectral width is 
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Another simple case in which equation (15) can be approximated in closed form is that in 
which the boundary flux is uniform, isotropic and has a Gaussian energy profile. The primary 



flux for this case is obtained by integrating equation ( 20 ) over all boundary points Xq and all 


directions H 0 . This yields the expression 
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By way of illustration, the case is considered in which the target is an aluminum half-space, 
z> 0, whose boundary is exposed to a uniform, isotropic flux of 56 Fe ions with a Gaussian 
energy profile having mean energy Eo = 1000 MeV/amu and energy width se= 10 MeV/amu. In 
view of the geometry of this problem, it is convenient to describe the unit direction vector, fi , 

by the spherical polar coordinates (1, 7? a ) that are related to its Cartesian components via the 


equation 
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where a is the polar angle and y is the inclination angle. Since there is axial symmetry about the 
z-direction, results are presented only for the case in which a = 0 deg. 

Figures 1-3 show the primary flux as a function of inclination angle, y , and energy at depths 
z = 0, 5, and 15 g/cm 2 of aluminum when a = 0 deg. The attenuation with angle and depth is 
clearly shown as is the increasing variation with angle at greater depths. 


THE PRODUCTION CROSS SECTIONS 

Before constructing the first order Green's function, it will be necessary to obtain an 
expression for the double differential production cross sections, a jk (Q&kiE,E k ) . The 
production cross section, <7j k (£l,il k ,E,E k ), is for fragments of mass Aj moving in the direction 



£2 with kinetic energy E that are produced when projectiles of mass A k moving in the direction 
£2* with kinetic energy Ek strike a target. In the analysis to follow, a tilde (~) will be used to 
distinguish between kinetic energy and total energy. Thus, if a particle has kinetic energy E, it's 
total energy will be denoted by E = E + m p , where m p = 938.272 MeV is the proton rest mass. 
The nuclear fragmentation cross section can be expressed in the form 

djt = a jk( E k)fjk( P)> (25) 

where a ]k {E k ) is the total cross section, and f ]k (p) the fragment momentum distribution. 

Experiments show that when fragments are produced, they may be projectile-like or target-like. 
This implies that the fragments have small momenta in their respective rest frames. It has been 
determined that the momentum distributions of these fragments can be well approximated by 
Gaussians in reference frames close to their respective rest frames. Here, it is assumed that the 
distributions are Gaussian in the frame where the average fragment momentum is zero. 
Quantities evaluated in this frame will be denoted with the symbol * to distinguish them from 
their counterparts in the laboratory frame. 

According to Townsend et al. [15], it may be assumed that the fragmentation is isotropic in 
the * frame, in which case the momentum distribution may be expressed as 
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where p* = p^ + p* , and the terms p x and pjf are the fragment's momentum transverse and 


parallel to the beam. The momentum width, cr, is given by Tripathi et al. [16] as 
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On transforming to the lab frame, while keeping in mind that /* transforms like —5- and 
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In equation (28), E = 7 ^-/^), pf = 7 l(P||-)M0’ , and(3 L = - y L ~ 2 . 
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The Lorentz factors, are given by 7 r = — for projectile-like fragments and r ) L = — — 
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for target-like fragments. The term E s = E s (E k ) is the lab frame energy shift which is 

determined as follows. 

According to Tripathi et al. [16], the projectile frame momentum shift is given by 
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and the corresponding energy shift by ^P s 2 + m p 2 — m p . On transforming to the lab frame it 


is found that 
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Let 0 = cos : (n • H *.) be the lab scattering angle and let (p,6,<p) be spherical polar 


coordinates with polar axis in the incident beam direction, Then, p± = psin0 and 

p\\ = p cos 0 , where p = ^ E 2 — m 2 . Since there is axial symmetry, 
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showing that the production cross section is given by 
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where 

E [E, 0] = 7 l [E - /3 L ,jE 2 - m 2 cos 0] . (33) 

Since the exponential appearing in equation (32) achieves its maximum when E = 7 m p and 

6 = 0 , a useful simplification can be obtained by expanding about this point in powers of the 
parameters e and p that are defined by e = E — 'ym p , p = 1 — cos 0 . By retaining terms up 

to second order, it can then be shown that 
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and 
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FIRST ORDER GREEN'S FUNCTION 


The next term in the Neumann Series that needs to be calculated is the first order Green's 
function, which represents the transport of first generation fragments. Since the Green's function 
boundary condition (10) is a special case of the broad Green's function boundary condition (17), 
we need only deal with the latter. It may be recalled that the first generation flux is given by the 

recurrence relation O 1 = (Q • L • S) • <X>° . Therefore, on replacing <t>° by the broad zero order 
Green's function and expanding the result, the broad first order Green's function is given by 
expression 
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where P \ " = x"- ft ! and Xj ' = x"— (pi M — P \ ')ft i is the point where the ray through x", in 


the direction ft i enters V. The residual energy is 



EM’-Pt'.E i) = + 


(38) 


and the mean energy at depth (pi pi ') of an m-type ion that entered the transport material 
with energy E 0 is 

4,(ft ft - (ft ft')]- (39) 


The term s h m (pi p 1 \E 0 ) is the broad spectrum energy width. 

The expression in equation (37) can be evaluated by numerical quadrature, but this is 
computationally expensive. Therefore, it is desirable to construct an analytical approximation. 
This can be done by making use of Taylor's theorem, the mean value theorem, and saddle point 
techniques as described in [10]. This yields the result 
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where p * is the root of the equation g jm (p ") = 0 , and the other terms in equation (40) are 
given by the following 

9jm(p ") = E m (p"~ P,E 0 ) - ^[i m (p"- p,£ 0 )] - Ej(p - p",E ), (41) 
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Ej(p — p ", E) = Rf^RjiE) + p- p"\, (44) 

E m (p"- P,E 0 ) = R m ~ 1 [R m (E o) - (/>"- p)]. (45) 

For the special case in which the boundary flux is uniform, isotropic and has a Gaussian 
energy profile, the first generation fragment flux is obtained by integrating equation (40) over all 
boundary points, x 0 , and all directions, Hq . This yields the result 
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Again, the results are illustrated for the aluminum half-space problem described above. 
Figures 4-6 show the secondary 16 0 fragment flux as a function of inclination angle, 7, and 
energy at depths z = 5, 10, and 15 g/cm of aluminum with an 1000 MeV/amu Fe projectile. In 
these figures, at a depth of z = 5 g/cm 2 , the number of fragments are increasing with depth until 
~ 80 deg where a large decrease with depth occurs as the stopping range of the particle is 
approached. Close to the boundary, there are no fragments produced, which is why there is still 
an increase in fragments at a relatively shallow depth. The other two figures show that with 
increasing depth, the number of fragments decrease as the stopping range is approached. 

As a verification of the approximations used, the first order fragment flux was computed by 
evaluating the integrals in equation (37) numerically. A comparison was then made between the 
numerical results and those obtained by analytical approximation. This comparison is shown in 
Figures 7-8 which exhibit the O fragment flux at depth z = 10 g/cm . Figure 7 shows the 
variation of the fragment flux with inclination angle, 7 , for several energy values and Figure 8 


the variation of the fragment flux with energy for several values of the inclination angle, 7 . Note 



that there is good agreement between the results, indicating that the approximations used are 
accurate. 

RESULTS 

In this work, analytical approximation techniques have been used to obtain closed form 
expressions for the first and second terms in the Neumann series solution of the three 
dimensional linear Boltzmann equation. Results were presented for the case in which a uniform, 
isotropic flux of 56 Fe ions with a Gaussian energy profile strikes the boundary of an aluminum 
half-space. A comparison between the closed form results and comparable results obtained by 
numerical quadrature shows that the approximations used are accurate. In addition, while the 
numerical results presented took several hours to generate, the analytical results were obtained in 
a few seconds showing the feasibility of obtaining a fast and accurate three dimensional ion 
transport code. 

The construction of the remaining terms of the Neumann series will be addressed in future 
work along with comparison to measured laboratory data. It is anticipated that an approximate 
closed form expression will be found for the second order Green's function and the series 
remainder obtained by a non-perturbative technique. 



NOMENCLATURE 


® = Flux column vector 
ft = Direction of propagation 
Q = Quadrature operator 
L = Linear transport operator 
S = Fragmentation or collision operator 
Pj(E) = Total survival probability 

Sj(E) = Energy lost per unit path length 
<jj{E) = Macroscopic absorption cross section 

E,E') = Double differential production cross section 
^.(x, ft', E') = Flux of k-type ions 

m p = Rest mass of a proton 7 1 = Fragment Lorentz factor 

/ 3 l = Fragment beta factor 

Eg = Fragment average energy downshift 


a = Fragment momentum width 
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FIGURE DESCRIPTIONS 


56 • 2 

Figure 1 . The Fe primary ion flux at depth z = 0 g/cm when a = 0 deg. 

Figure 2. The Fe primary ion flux at depth z = 5 g/cm when a = 0 deg. 

Figure 3. The 56 Fe primary ion flux at depth z = 15 g/cm 2 when a = 0 deg. 

Figure 4. The O fragment flux at depth z = 5 g/cm when a = 0 deg. 

Figure 5. The 16 0 fragment flux at depth z = 10 g/cm 2 when a = 0 deg. 

Figure 6. The O fragment flux at depth z = 1 5 g/cm when a = 0 deg. 

Figure 7. A comparison between the numerical and approximate 16 0 fragment fluxes at depth 
z = 10 g/cm for several energies (MeV/amu). 

Figure 8. A comparison between the numerical and approximate 16 0 fragment fluxes at depth 
z = 10 g/cm for several inclination angles (deg). 
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